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SUMMARY 


A method for analyzing the nonadiabatic viscous flow through 
turbomachine rotors is presented. The field analysis is based 
upon the numerical integration of the full incompressible stream 
function-vorticity form of the Navier-Stokes equations, together 
with the energy equation, over the rotor blade-to-blade stream 
channels. The numerical code used to solve the governing equations, 
employs a nonorthogonal boundary-fitted coordinate system that 
suits the most complicated blade geometries. The effects of 
turbulence are modelled with the two equations which will be 
reported in the second volume of this contract. A numerical 
scheme is used to carry out the necessary integration of the 
elliptic governing equations. 

The method of analysis is applied to various types of turbo- 
machine rotors. First, the flow characteristics within the rotor 
of a radial inflow turbine are investigated over a wide range 
of operating conditions. The calculated results are successfully 
compared to existing experimental data. Second, the flew in a 
radial compressor is analyzed in order to study the behavior of 
viscous flow in diffusing cascades. The results are compared 
qualitatively to known experimental trends . The solution obtained 
provides a great insight into the flow phenomena that takes place 
in this type of turbomachines. Comparison with nonviscous flow 
solutions tend to justify strongly the inadequacy of using these 
solutions with standard boundary layer techniques to obtain 
viscous flow details within turbomachine rotors. It is concluded 
that the method of analysis is quite general and gives a good 
representation of the actual flow behavior within turbomachine 
passages . 

The computer used in this work is an AMDAHL 470. The flow 
domain has been divided into 30 step sizes in rj direction and 
40 in the £ direction. Typical CPU time was 120 seconds. 


1 


INTRODUCTION 


For many years it has been recognized that the flow in turbo- 
machines is characterized by the presence of three-dimensional r 
viscous, and compressible effects occurring in a complex 
geometrical configuration. The basic understanding of the physical 
phenomena in this flow requires, therefore, an analysis involving 
the solution of the unsteady three-dimensional, compressible viscous 
flow equations within the rotating and stationary blades 
comprising the machine. Such an analysis is clearly a formidable 
task. The complexity of the entire problem requires some kind of 
simplification for one or more of the important factors affecting 
the flow. These simplifications must cover the essential physical 
process with sufficient quantitative accuracy and still permit a 
clear and rational calculation of different flow processes. Most 
efforts have been concentrated frequently on the solution of the 
steady and inviscid version of the flow governing equations. 

These solutions have been marked by increased versatility in the 
ability to deal with subsonic [1, 2, 3 and 4] as well as transonic 
flow regimes [5 and 6] . In order to give a more accurate repre- 
sentation of the actual flow processes, approaches [7 and 8] have 
been devised to account for viscous effects. Most of these 
approaches are based on the assumption that a two layer model is 
representative, i.e. an inviscid flow solution which interacts with 
an end wall boundary layer solution. Important contributions to 
viscous flow analyses in turbomachinery have been made more 
recently' in References [9, 10 and 11]. The important features of 
these are the attempts to solve the parabolized version of the 
complete three-dimensional viscous flow equations with special 
techniques . 

Although the above remarks are not intended to be a complete 
survey of all the available methods, it is evident uhat the inviscid 
analysis is useful for providing a considerable insight into the 
character of the flow. However, the neglect or viscous effects is 
a serious shortcoming if detailed quantitative information is 
desired to calculate viscous losses or heat transfer. The approach 
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used to account for viscous losses effects by using the boundary 
layer technique proves to have a number of drawbacks. First, the 
correct means for matching boundary layer and inviscid solution 
has not been established if the inviscid flow is rotational, such 
as in the case of a curved or a rotating passage. Second, most 
of the existing interacting boundary layer analyses are not capable 
of handling strong interaction mechanisms of the types present in 
turbomachine rotors. The parabolic flow approximations, on the 
other hand, neglect completely the downstream influences. Con- 
sequently, important effects such as surface curvature, downstream 
blockage and reversed flow regions are totally ignored in these 
type of approximations. To circumvent this deficiency, a procedure 
based on the solution of the full elliptic Navier Stokes equations 
is required. Unfortunately, such direct procedures have defied 
accurate numerical solutions due to the limitations imposed on 
the core size and speed of present computers. Moreover, the lack 
of powerful numerical schemes capable of achieving a rapid con- 
vergence for three-dimensional elliptic equation renders the 
solution to be costly. 

In the present study, an attempt is made to demonstrate the 
feasibility of obtaining viscous flow details within turbomachine 
passages by appropriately combining several blade-to-blade viscous 
flow solutions. Each of these solutions is obtained through the 
numerical integration of the full Navier-Stokes equations over 
a predetermined computational surface that extends between the 
blades. The set of computational surfaces required for the analysis 
are themselves generated from the solution of the nonviscous 
version of the Navier-Stokes equations, as suggested by Wu [1] . 
Because of the constraints implied by the use of these computa- 
tional surfaces the resulting viscous flow details are regarded 
as a quasi- three-dimensional description of the flow field. The 
use of a non-viscous flow solution to generate a viscous solution 
in the manner outlined is anticipated to overcome some of the 
drawbacks of the method discussed earlier. For example, the 
inclusion of the full Navier-Stokes equations in the solution 
procedure makes it possible to account, in an effective manner. 
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for the interaction between the viscous and inviscid regions. 
Moreover, preserving the ellipticity of the problem by working 
with the full equations, besides offering an accurate repre- 
sentation of the flow field, allows for the recognition of down- 
stream influences. 

The study will be presented in two parts. The first part 
which is reported in the present volume deals with the general 
f ormulation of the viscous flow equations , and presents the 
results for laminar flow cases. Turbulent flow cases considered 
in this investigation will be reported in the second volume. 

The present volume consists of four sections described 
as follows. The equations that govern the flow of viscous 
fluid within turbomachine passages are presented along with a 
closure model that accounts for turbulence effects in section 1. 

A rigorous discussion regarding the accurate representation of 
the different boundary conditions is also given with special 
emphasis placed on the determination of the downstream boundary 
conditions that are required to preserve the ellipticity of the 
problem. A transformation of the flow equations using a non- 
orthogonal boundary fitted coordinate system which is numerically 
generated using Thompson Code [12] is presented in section 2. 

The overall effect of this transformation is to produce a domain 
in which the arbitrary blades shapes of the turbomachine become 
straight and parallel. The details of the numerical scheme used 
to integrate the flow governing equations as well as the pro- 
cedures employed to handle the nonlinearities in these equations 
are given in section 3. Finally, the results of the application 
of the present method of analysis to various types of turbomachine 
rotors are presented in section 4 . 
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1. FLOW ANALYSIS 


This section presents the detailed development of a method 
for analyzing the viscous nonadiabatic flow of gas through 
turbomachine rotors. The partial differential equations that 
govern the flow behavior within the machine passages are presented 
first. A transformation of the general equations from the 
three dimensional form to several particular two dimensional 
forms, on predetermined stream surfaces, is then outlined. 

The resulting flow equations are further expressed in the con- 
servation law form using the vorticity- stream function formulation. 
This is followed by a discussion of the necessary boundary condi- 
tions that lead to a unique solution to the problem. 

Fundamental Aerothermodynamic Relations 

The three dimensional viscous, compressible flow within 
turbomachine rotors is governed by the following set of laws. 

Conservation of Mass ; 

+ v* (pW) = 0 (1) 

Conservation of Momentum : 

Newton’s second law of motion when combined with Stokes 
hypothesis can be written as 

- 2 

p(|~ + {W*V)W + 2flxW - ~ vr 2 } = - Vp - Vx[y(V x W) ] 

+ i V(pVW) (2) 

where fi is the rotor angular velocity. 

Conservation of Energy : 

In the absence of heat sources or sinks, the first law of 
thermodynamics for a fluid with a thermal conductivity, K, can 
be written as: 

p(f| + W*Vh) = ||: + (W-V)p + D + V * (KVT) (3) 
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In the previous equations, p, p, T, and h denote the static 
pressure, density, temperature and enthalpy, respectively. 

While, u represents the kinematic viscosity coefficient and 
D is the viscous dissipation function. W is the relative velocity 
vector at any point, whose location is defined by the relative 
position vector ?, in the rotating frame of reference. The 
relations among the flow state variables are those of the ideal 
gas 


p “ 

pRT 

(4) 

dh : 

= C dT 

(5) 


P 


The mass-averaged variables principle is utilized in the above 
equations in order to describe the behavior of a turbulent flow. 
According to this principle, the components of the velocity 
vector W, and the temperature T, are expressed on a mass-averaged 
basis, while the pressure p and tne density p are expressed by 
their mean values only. For any flow variable O (x-^ , , t) , 
its mass-averaged value, q, is given by: 

q(x 1 ,x 2 ,x 3 ) * Q(x 1 ,x 2 ,x 3 , t) p* (x 1 ,x 2 ,x 3 ,t)/p 

where p* (x^ ,x 2 ,x 3 , t) is the instantaneous value of the density 
at any point whose location is defined by the coordinate (x^,x 2 ,x 3 ) . 
The over bar in the above equation represents the conventional 
time average, thus p* s p. An effective viscosity u is also 
used instead of the kinematic viscosity u . This effective 
viscosity is assumed to describe the effects of Reynolds stress. 

In the present work, a tx*o equation turbulence model will be 
employed to calculate The complete details of the model 

and its implication is outlined in Volume 2 of this study. 

Stream Surface Equations 

An approach will be taken to reduce the spatial dimensions 
of the problem from the general three dimensional form to several 
particular two dimensional forms. The objective is to obtain a 
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solution to the flow governed by the equations (, 1) to (5) through 
an appropriate combination of two dimensional flow solutions. 

The reduction in spatial dimension is achieved through the 
consideration of the blade-to-blade stream surface concept. 

The blade-to-blade stream surface, S^, may be described by the 
annulus that would extend from the pressure surface of a blade 
to the suction surface of the next blade, as shown in Figure 1. 
This annulus is characterized by the variation of its filament 
thickness, b, and the radius, r. For the purpose of the present 
discussion, the stream surface will be considered to represent 
the mean geometric properties of the annulus. It is possible 
to trace out the shapes and the filament thicknesses for a finite 
number of these blade-to-blade stream surfaces (annuli) in any 
turbomachine passage, using a meridional flow analysis as that 
indicated in reference [1] . On these stream surfaces the flow 
equations (1) , (2) and (3) are transformed to several two 

dimensional mathematical expressions. This is possible since 
each stream surface provides relations among the coordinates, 
such that the variation of flow properties over each surface may 
be described in terms of two space variables only. If a solution 
to the resulting two dimensional equations is obtained for each 
stream surface, then the flow properties throughout the three 
dimensional field may be readily evaluated. 

Considering the flow annulus shown in Figure 1. The 
curvilinear distance along the intersection of the mid-surface 
of the annulus with a meridional plane is denoted by m. The 
distance normal to the mid-surface is represented by n. The 
circumferential coordinate <j> is considered positive in the. 
counterclockwise direction when viewed down the positive z axis. 
The thickness of the annulus, b, is assumed to be small compared 
to the radius r. Hence, the n component of the velocity vector 
and all variations in the n direction are neglected. The trans- 
formation of the flow governing equations (1) and (2) to the 
stream surface (annulus) coordinate system (m-<j>-n) is outlined 
in Appendix A with the following results . 
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Continuity : 


If + Is < brp V + I? tbp V - 0 


(6) 


Meridional Momentum: . 


3W 


3 W W, 3 W w 2 
m . 4 m Wo _ • 


P + p(W m 3if + r 3T"— 5ina “ fl r sina “»V ina) 

3£ 2 3 aina + 1 8H *n 

= " 3^ "* 3 3^ [tJ e { 3 IT T “ Sina + r 34“ } 1 

- - 3W a 3W, W, , 8W 

+ I <k <2u e r 35r )+ sy'Vair-s 1 sina + ? rr 11 


y a , 3 W, w ' 

- 2 — [— 7 -r^- + — sina] sina 
r r 34 r 


( 7 ) 


Tangential Momentum : 

aw, aw. 

pr rr^- + p (rW 3 -^ 
^ 3 t m 3m 


3W, 

+ W — * 


■ - 5 ^— + W W, sina + 2flr W sina) 

4 34 m 4 ) m 


, * 3W W , aw. 

= - ? si ft* (rrr-+ — sina + - ] 


34 3 34 Lt e'3m r 

_ 3W. 3W « -, 3W W 

+ to tlJ e !r V iba+ 3ir n ' f 3? I2u e'F 5F” + — sina)I 


r 34 


3W , W 


+ “ e sin “ 'sr 


4 * , 1 3 t.t T 

sina + — r-r W ] 
r r 34 m 


(8) 


where W^, are the components of the mean relative velocity 
vector W in the meridional and tangential directions, respectively, 
and a is the angle between the mid-line of the annulus, m, and 
the axis of rotation, z. 

The system of the differential equations ( 6 ) , (7) and ( 8 ) , 

can be written alternatively in terms of a stream function 4 
and a vorticity w. In the present study, the stream function- 
vorticity formulation is considered for the purpose of reducing 
the difficulties of the coupling and nonlinearities associated 


with the presence of pressure in the flow equations (7) and (8) . 

In addition, this formulation offers the possibility to express 
the governing equation in their conservative law form. Con- 
sequently, no excessive accumulation of errors in the fluxes of 
the conserved quantities will result in the finite difference 
approximation of the governing equations. Roache [131 
illustrated this by showing how the Gauss divergence theorem is 
satisfied for the finite difference equations. 

For steady state flow, the continuity equation (6) is satisfied 
by introducing a stream function , which is related to the flow 
velocity components by: 


w = m i_ 
m b pr 3(j) 


and 



M 1 3il> 
b p 3m 


(9) 


where M is the mass flow rate passing through the annulus, S^, 
of Figure 1. If the pressure terms are eliminated from equations 
(7) and (8) by cross differentiation and the mean vorticity 
variable m is introduced, one obtains the vorticity transport 
equation: 


<2 II «) - It Cir »> - 1= tr h; 1 


.M 3tji 


3m 'b 34> 


3 'b 3 m 


3m 


3m 


+ h r r h (y e u) 1 + G 1 " 0 


where W“ = t W? 

m <p 


G l = 2!U !h ( E sin “ H )_ I? ( l sin “ 


. 3p 3WV2 
3m 3(f) 

and u is defined by: 

“ = r l fm (r V 


3p Spl + r slno |j 

3 3m 3<f> 


f? « V 1 


( 10 ) 

(ID 


( 12 ) 


(13) 
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T 


When the velocities W m and W^ in equation (13) are expressed in 
terms of the stream function variable, as defined in equation (9) 
equation (13) would reduce to the stream function equation: 


1 . 3 . Mr 3<|n 3__ , M 3 tJj , •, 

r L 3m l bp 3m J 3$ ^bpr 3<jr J 


(14) 


Energy Equation 

For a turbomachine rotor it is convenient to express the 
energy equation in terms of the total enthalpy (H) of the gas, 
besides its velocity components. The total enthalpy for 
turbulent flow is expressed as follows: 


H 


Cp T + 



+ QW.r 
<p 


, n 2 r 2 
+ -s + E 


(15) 


where E is the kinetic energy of turbulence. Thus, the energy 
equation (3) , when transformed to the stream annulus coordinate 
system, as given in Appendix A, can be written as: 


I l k (H If* • f* !H 1 


3,^e 3H, 1 3, He 3H, 

3m l Pr 3nr r 3<{> 4 Pr 3<j> 


+ 4={u„r S§n 


3m 


Pr 3m 


k S CE Pr' 3m 


3 ri- 3W 2 /2 ,1 1_. 3E. , 

34> l r l Pr 3<f> " ^S CE Pr ; 3<p Ji 


# f? <“e“> - Dr + G 2 = 0 


"♦ IS 


(16) 


where Pr is the turbulent Prandtl number and S__ is the turbulent 

Cllj 

Schmidt number for the kinetic energy of turbulence, E. The 
source term in the above equation represents the generation 
or decay in energy, due to the effect of rotation. It is given 
by: 


g 2 - - n I ( fs'<v +r2al H ] - l^r+A) Ml) 


4* < 7 {^ — -r iL— (yf >- 4 . \ ] -l i 3_..r e 3 /r,j rr \ ] i (17a) 

41 3m l pr “ 3m 1 6 r 2 ‘ J ~ ~ tw^r+ — ; J } 


r 3<J) Pr 34) v v 


The dissipation function, D, can be written as 

aw 2 - aw, w _ 

D= 2w e {< aif) + <lsi i + r lsi ““) > 
aw, 7 aw w, „ 

+ W e { a sr ' rW~ ~ slna> 


(17b) 


The properties of flow passing through the stream annulus S^, 
are completely defined by equations (4) , (9) , (10) , (14) and (16) 
together with the known variations of y g , Pr and the given boundary 
conditions. The effective viscosity, u , is calculated from a 
two equation model, one expressing the development of the turbulent 
kinetic energy, E, and the other its dissipation rate, e. These 
equations may be expressed in terms of the stream annulus system 
of coordinates (m,<j),n) as follows^ 


Turbulent Kinetic Energy Equation : 

a 


« [1_( E 1 £* 
b 3m * 5 aV 


3$ 


3<j>' 


/r* ftiih l _ ° 

3(f) 3m' J 3m 


b CE 


5 CE 


3E 

r3<f> 


} - rD + per = 0 


r Mi 
3m ; 


(13) 


Dissipation Rate Equation; 


M 

b 



3 


r — ( Ue _ r i£) 

3m V S_ 3m ; 

L, c. 


5- 

3c l S 

Ce 


I i£). c 

r 3<fr 


£ 

E 


rD + C. 


£_ 

E 


pr = 0 


(19) 
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where S Cg is the turbulent Schmidt number for the dissipation of 
kinetic energy of turbulence, e. The values of the constants 


C l' C 2' 
Table 1 . 


C D r &CE an ^ S Cs fc ^ e tur ^ u l ence model are given in 


TABLE 1. 

VALUES OF 
k-£ MODEL 

THE EMPIRICAL 
OF TURBULENCE 

CONSTANTS 

FOR THE 

* 


c„ 

s„„ 

S„ 

D 

JL 

2 

CE 

C£ 

0,09 

1.42 

1.92 

1.0 

1.3 


The expession for the effective viscosity y is given by 

s 


v 


e 




+ C 


D 



( 20 ) 


where the laminar viscosity, y^, is considered in the present 
study to be uniform and known. 

In general, equations (4) through (20) are valid for any turbo- 
machine geometry or any number of stream annuli except for the two 
stream annuli Sq Q and Sq^. shown in Figure 2, which contain the 
hub and shroud contours. This exception may be attributed to the 
existence of a large variation in flow properties along the nor- 
mal, n, to these two annuli resulting from the presence of the 
solid boundaries. The determination of the flow properties within 
these stream annuli constitute a study by itself and is not 
intended for inclusion in the present work. 

The solution of the above system of equations within the 
turbomachine passages is carried out numerically. One can observe 
that equations (10) , (14) , (16) , (18) and (19) , constitute a 
system of coupled elliptic partial differential equations, involving 
second order derivatives of ty, to, H, E, and e which are the dependent 
variables. From the nature of the problem, none of the terms are 
negligible in the governing equations. The convective terms 
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introduce nonlinearity and also instability if. the proper dif- 
ferences are not taken into account. Once a solution for these 
variables has been obtained , the~velocity distribution can be 
determined from equation (9) . The pressure distribution can 
then be evaluated from either equations (7) or (8) . 

Description of the Computational Domain in the Physical Space 

In order to solve the elliptic equations by the usual numeri- 
cal methods, it is necessary to define a region in the physical 
domain with boundary conditions specified for the different 
dependent variables on all boundaries. The flow region of interest, 
as shown in. Figure 3, contains the blade row and segments of the 
stream surface, S^, extending upstream and downstream of the row. 

Due to the circumferential periodicity in the turbomachine 
passages, the selected domain need to encompass only a fraction 
of the flow annulus containing a single blade to blade passage. 

The shape and location' of the periodic boundaries (AB, NM, IH 
and FG) may be defined arbitrarily as long as their spacing 
corresponds to the blade pitch. The upstream and downstream 
boundaries (AN, GH) are located sufficiently far from the blade 
so that the tangential variation in flow properties along them 
is negligible. The flow properties are consequently considered 
to be uniform along the boundary AN and GH. 

Boundary Conditions 

In specifying the boundary conditions, two flow cases are 
investigated. Preliminarily, only the case of the laminar flow 
is considered in this report. The turbulent flow case will be 
reported in a second report. Accordingly, in the following 
specification of boundary conditions, no assignment for the boundary 
values of E and e in equations (18) and (19) is needed. Moreover, 
the flow properties within the turbomachine channels will be 
completely defined through the simultaneous solution of equations 
(10), (14) and (16). 
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a. The Upstream Boundary AN : 

It is a common practice in turbomachine flow calculations that 

the magnitude and direction of the flow velocities , the total 
* 

temperature, and the total pressure or density are defined at 
the turbomachine inlet. Therefore, along the boundary AN, the 
values of ij>, w, H or their derivatives can be evaluated using 
the defined flow properties. The known magnitude of the inlet 
relative velocity and its direction, as shown in Figure 3, 
specify the values of 3i|>/3m and according to the following 

relations: 


M = 

3m 


bp 

M 


"♦ “ 


r 3$ tanS inlet 


_ bpr w _ _1 

9<f> m m 2lf / z 


( 21 ) 


where Z is the number of blades. Since the inlet stream of gas 
is considered to be uniform, the absolute value of the vorticity w 
has to be zero along the boundary AN. In a rotating frame of 
reference, the relative value of to is given by the following 
expression: 

to =« - [20 sinct] in iet ( 22 ) 

The value for the total enthalpy, H, can be defined using the 
specified flow properties. 

b. The Periodic Boundaries AB, NM and FG, IH t 

The periodicity condition requires that the direction and 
magnitude of the flow velocity as well as other fluid properties 
be equal at every two corresponding points along AB and MN. 
Similarly, the same conditions should apply at every two corre- 
sponding points along FG and IH. In terms of the present 
dependent variables, the periodicity requirements are satisfied 
through the following conditions. First by equating w, 3to/3d>, . 


14 


3ij>/3<i>, H and 3H/3<j>; values at each two corresponding points. 
Second, ensuring that the ijj values differ by unity between the 
corresponding points. 

c. The Blades Surfaces Boundaries Ml and BF : 

For the laminar flow case, two boundary conditions over 
the blade surfaces are usually specified. These are the non-slip 
condition and the impermeability of the surface in the case of 
blades with no injection. The non-slip condition requires that 



where N is the normal to the blade surface. On the other hand 
the impermeability condition requires the component of velocity 
in the direction normal to the blade surface to vanish. Therefore, 
the blade surfaces are treated as streamlines with the \ p values 
specified as zero on the MI surface and unity on the BF surface. 

On either MI or BF surfaces, one has therefore, two boundary 
conditions for if) but none for w. It is a well accepted fact in 
computational fluid mechanics to rely on a modified evaluation 
of equation (14) to determine the boundary condition for cu . 

The modification is introduced in an attempt to insure that 
equation (23) holds, that is, to satisfy the no slip condition. 

This approach is utilized in the current study to determine the 
value of the vorticity, w, over the blade surfaces. The details 
of the procedures used will be presented in the next section. 

In regard to the thermal boundary conditions either the 
blade surface temperature is known or the normal derivative 3T/3 ti 
is specified as zero for the adiabatic wall conditions. In 
either case, equation (15) is used to determine the value of H 
or its derivative along the blades surfaces. 

d. The Downstream Boundary GH ; 

A few basic problems arise in the specification of the 
boundary functions for the dependent variables along GH. The 
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first concerns the bheavior of the dependent variables H and as. 
The nature of the problem is that, physically, these are known 
only if the boundary GH is located at an arbitrarily large 
distance from the blade surface. In this case the value of H 
and tu is that of the corresponding surroundings. The placement 
of GH at exceedingly large distances from the blade boundary is 
quite obviously not possible for numerical considerations. 
Therefore, one has to employ some auxiliary conditions, usually 
obtained by experience, to define H and os implicitly. The 
conditions of zero gradients in the meridional direction, m, 
is employed in the current work. 



For the vorticity, to, the absolute value is taken to be zero, 
hence , 


os = - [2£2sinct] ^ (24b) 

More important than the specification of the remote boundary 
functions of to and H along GH is the determination of the values 
at the same boundary. The downstream flow velocities, which may 
be used to determine the stream function derivatives along GH, 
and that guarantee a unique solution to the problem are not 
known in general apriori. Therefore, one has to introduce a 
supplementary condition, generally resulting from physical 
intuition, to define the stream function derivatives. Investigators 
working in the inviscid flow area dealt with this problem by 
using an iteration procedure, through which the Kutta condition 
for tangency of the flow at the blade trailing edge is satisfied. 
This is equivalent to specifying a unique solution to the problem. 
Unfortunately, the Kutta condition cannot be applied realis- 
tically to solve the present flow problem due to the viscosity 
effects. The conservation of angular momentum principle [14] 
is employed as an alternative supplementary condition that results 
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in the required unique solution. The iterative procedure used 
to implement this condition is as follows. 

Estimated exit flow angles, B a ^ on 9 GH are used to 

specify the values of stream function derivatives in the m 
direction through the following relation: 


3 £ _ _1 

3m 2 tt/Z 


tan 6 


exit 


‘exit 


(25) 


The flow field equations are then solved for the boundary function 
of given by equation (25) to obtain the velocity and the 
pressure distribution throughout the stream annulus of Pig. 1. 
An evaluation of the torque developed by the annulus is obtained 
through the integration of the pressure and shear forces acting 
on the blade surfaces. The change in the angular momentum 
between the known inlet and the estimated exit flow conditions 
is determined. If the value of the predicted torque was not 
equal to the rate of change of the angular momentum, then the 
direction of the exit flow velocity is altered. The whole pro- 
cedure is repeated until a satisfactory result is obtained. 

The solution of equations (10) , (14) and (16) , subjected 

to the above boundary conditions within the blade rows and in 
the near field, are carried out numerically. In this connection, 
it is necessary to reduce the complexity of handling the finite 
difference representation of the governing equations and the 
associated boundary condition near the curved boundaries of the 
blade surfaces. This is best accomplished by introducing a 
coordinate transformation from the (m-4>-n) system of coordinate 
to a contracted body-fitted coordinate system. The overall effect 
of this transformation is to produce a square field in which the 
arbitrary blade shapes become straight and parallel. The 
development of such a coordinate transformation is presented in 
the next section. 
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2. BOUNDARY-PITTED COORDINATE SYSTEM 


In all fields cone- rned with the numerical solution of partial 
differential equations , the physical region in which the solution 
is desired is overlayed with a grid. In constructing a grid over 
the blade-to-blade configuration of Fig. 4, the points on the 
blade surfaces do not generally correspond to grid points. 
Interpolation must therefore be used to determine the function 
values immediately adjacent to the boundary points, for the given 
boundary conditions. Moreover, if Neumann type boundary conditions 
are present, interpolation is also required to determine the 
boundary values themselves. Interpolation between grid points 
not coincident with the boundaries is particularly inaccurate 
in the case of field equations that produce large gradients 
in the vicinity of the boundaries [13] . This inaccuracy in 
representing the boundary conditions is known to impair the 
success of any numerical scheme in achieving an accurate con- 
vergent solution [15] . It is therefore desirable to use a 
coordinate system such that the problem boundaries lie along the 
coordinate directions. Such coordinate is commonly defined as 
a boundary-fitted coordinate system. 

In this chapter, the available methods for developing 
boundary fitted coordinate systems for general shaped bodies 
are briefly discussed in the first section. The procedure used 
to transform the physical domain of Pig. 3, to a unit square 
using the boundary- fitted coordinate system is then outlined in 
section ii. The flow governing equations with their associated 
boundary conditions in the transformed domain are presented in 
the third section. The line integral method employed to obtain 
the pressure distribution within the flow field is covered in 
section iv. 

Basic Transformation Methods 

The importance of generating a boundary fitted coordinate 
system in viscous flow problems is evidenced by the fact that 

i 

the only successful Navier-Stokes solutions to date have been for 
those bodies for which such coordinates is available. For 
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simple geometric shapes, it is usually possible to employ simple 
algebraic transformation to place one of the coordinates along 
the boundary surfaces. However, for more complex boundaries 
such as those involved in turbomachinery applications, it is 
extremely difficult if not impossible to use an analytical 
treatment to generate a boundary-fitted coordinate system. In 
these cases the boundary fitted coordinate system is generated 
numerically . 

In reference [16] , Stanitz takes the boundary-fitted coor- 
dinate to be the stream lines and equipotential lines that 
result from the solution of Laplace equation for the ideal two 
dimensional flow over the area of interest. Although this 
approach is straightforward and simple, it is strictly limited 
to two dimensions and is not particularly flexible in terms of 
coordinate spacing. A much more general method is to generate 
the boundary fitted coordinates by solving a pair of Poisson 
elliptic partial differential equations with Dirichlet boundary 
conditions. The boundary conditions specify one of the coor- 
dinates to be constant on each of the physical boundaries. A 
chosen distribution of the other coordinate is specified around 
the boundary contours. This procedure causes some coordinate 
lines to be coincident with each boundary of the physical domain 
regardless of its shape. The basic concept of such procedure 
has been employed in varied form by several investigators [18 , 
19,20]. Thompson [12] has extended this technique recently 
to be applied to multiconnected regions with any number of 
arbitrary shaped bodies. His method offers also the advantage 
of a provision for controlling the spacing of coordinate lines 
near any designated surface. These factors led to choose 
Thompson’s approach to generate a boundary fitted coordinate 
system for the blade-to-blade domain of Fig. 3. 
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ii. Mathematical Formulation 


Two transformations are employed in the present study to 
implement the generation of the boundary- fitted coordinate 
system for the blade- to-blade domain of Fig. 3. The first one 
is obtained by defining a stretched meridional coordinate , x, 
given by: 



This coordinate stretching maps the physical space of Fig. 5a 
onto the domain shown in Fig. 5b with (x, (*>) as coordinate 
system. 

The second transformation [12] generates the boundary 
fitted coordinates £ and n through the numerical solution of 
the following equations for x(£/n) and <p(£,ri) • 


J i!| . 2B a 2 * 
3£ 2 


353 n + Y ri " Q(5 ' n) 

oT] 


. J. a 2 * - n/r 5 

3£ 2 9 ^ 3n 3n 2 


Subjected to the following boundary conditions: 


x = q l 

(?/n x ) 

on 

ABFG 

0 = g x 

(^n 1 ) 

on 

ABFG 

x = q 2 

(S,n 2 ) 

on 

NMXH 

4> = g 2 

(S,n 2 ) 

on 

NMIH 


(27) 


(28) 


where 
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3x 

3n 


2 


) - '!£> 


a _ 0* 

P 5 V, 


3x 3x t 34i 3$ 
3n 35 35 3 ti 


(121) 2 + ( li) 2 

^35 J l 3? J 


( 29 ) 


The functions Q and P in the above equations are appropriately 
chosen to provide control over the spacing of the coordinate 
lines in the field. On the other hand, the functions q^, q^r 
g^ and g^ are usually specified by the known shape of the 
contours ABFG and NMIH in the (x, <j> ) domain. 

Mapping the region of interest in Fig. 5b in terms of the 
new boundary- fitted system of coordinates (5,n) yields a fixed 
square field in the final transformed domain as shown in 
Fig. 5c. Note that the blade surfaces in this transformed 
domain become straight and parallel. With this procedure 
the numerical solution of the flow equations developed in 
Section 1, is carried out on the fixed square field of Fig. 5c, 
using a uniform grid with no interpolation required regardless 
of the blades shape in the physical space. The transformation 
of the governing flow equations from the physical space to the 
transformed domain is outlined in the following section. 


iii. Transformed Governing Equations 

The general transformation from the physical space (m,<f>) 
to the transformed field of (£,n) is given by equation (26) 
together with the following vector function; 


5 


5(x,<|>) 

n 


0 (x,<j>) 


( 30 ) 
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This vector function as well as its inverse transformation are 
defined Once the (5,ri) system has been obtained. Partial 
derivatives of any scalar function, f, are transformed utilizing 
the chain rule as follows: 


i in 


at i / n if a$ sf, /T 

3m r l an 3^ “ 35 3n ;/ 


3f _ ,3x 3x 3 f . 

3<j> l 35 3n “ 3n 35 J/ 


{ 31 } 


where J is the Jacobian of transformation, given by: 


t — 3x 3 £ _ 3x 3<j> 

“35 3ti 3n 35 


Higher order derivatives as well as the derivatives normal to 
the different boundaries are presented in details in Appendix B . 
Using the expressions providing by the relations (31) and (B4) 
through (BIO) , the governing flow equations can be written as 
follows in terms of the new variables (£,n) : 


Vorticity Transport Equation 

3 ,M 3ifp . 3 ,M 3^i . 1 r . 3 2 3 2 , . , 

35 b 3n 3n ( b 35 5 J [6 a? 2 (y e u>) . 20 353t {y e w) 


' l: 


+ y 7~2 ^ + 5 In ( V“ > + T If + G 3 " 0 

3n (32) 

where the source term, is given by: 


G = i£ awV 2 - is. wr/2 + 2 o r a_ ,m sina 

3 5 3r| ari ac- + *«tar Vk sinct 


an 35 


35 'fa 


3£ 

3n 


. 3 ,M . 3ii>. . 


+ nV sino «S? H H H> 


(32a) 
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And the coordinate transformation parameters a and r are given 
by the following expressions: 


a = 


T = 


r3ji>.r s 

{ n [0 




- 3 2 x 






- 28 

3 2 X 
3 53 $ 

+ Y 

3 2 X , 

an 2 ” 

r- 28 

3 2 cb 

_L v 

a 2 *, 

an 2 

353n 



9x r 8 2 4) 
^ 8? 


-28 


2 

3<j> r 3 x , fl 
- 2S 


3 2 <j> 

353n 

3 2 x 

353$ 


+ Y 


3 2 d> 

3n 2 

3 2 x 


3n 

(32b) 


It is not difficult to show that the right hand sides of the last 
two equations reduce to zero if a uniform mesh size is used 
along the $ direction. In this case 7 the coordinate system 
will be referred to as a non- contracting coordinate system# 


Stream Function Equation 


5 (£L_ _ c S. ( L - q « / -- * y \ t -v ° f li "fi 

3? ^Sp d£ } B 3C bp 3n } ° 3p ( bp 3£' + Y 3rT { b? 3 t^ 


a*. 


at, 


3 m at, 


, rrr 3^, _ 2 _? 

* bp 3n + T 3c } “ r J 


(33) 


Energy Equation 


II [ 3 /jj _ 5 m \ i rx 3 . 3H , g 3 ,^e 3H % 

b l 3C tH 3ri 3n CH Tf }i j [5 aT { p? 3l } ~ s H ( p? 3n } 


an' 3n 


35' 


B I- v L ^ 3H . ^e , 3H 3H, . 

3n ( Pr 3£ I + Y 3n f p r 3 t^ + Pr ^ 3n + aT ^ 


§ + G 4 - 0 


(34) 


where 


■]}/J 
■3 >/J 
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If 
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+ "a {( k 3 C 


aw 2 aw. 
o © in , . . 3 x p 

35 an J l 35 9n 
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(34a) 


' (34b) 


( 34c) 


Turbulent: Kinetic Energy Equation 


M r9_ fp £&_ 3, p 3£ n „ If, 3 

b L 35 CE 3^ 3n (E 35 } 1 


f s J3E, o 3 , e 3E* 

a? l s CE 35^ aT's^T W 
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V, 


3 n v s CE 35 


u. 


3 f e 3E. 3, H e 3E, . -e 

^ c ^r) + T a-n l 


3 "' S CE 3n ) + S CE 


3E , 3E,, 
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D 2 

- j + Jp r e = 0 


( 35 ) 


Dissipation Hate Equation 


M rl_ (e 3ij_ i_ fe _ lr, 3_ , U e 3e. fl 3 /e 3e, 

b l 35 { 3n 3n ( 35 )] J 15 35 { S^ 35 )-8 35 ( S^ SrT 5 


Ce 


Ce 


8 de \ 1 -.- 3 / _e 3 e u y e 

3n ( s CE 3C> + y 3n ( s ce H ) + 


( rc ^ £ t 9 £ \ 

<0 3^ + T 3?> 


- e D - e 2 

C 1 E J * C 2 P E J r “ 0 


(36) 


The above system of equations (32) - (36) are somewhat more com- 
plicated by the extra terms added by the transformation. The 
disadvantage of having these terms, however, is far outweighed 
by the computational advantages of the simple square flow region. 

In general, one can demonstrate that the transformed flow equations 
are still elliptic in spite of the appearance of the cross deri- 
vative terms ^ . The numerical solution of these equations 

using a uniform rectangular gird in the (5,n) domain, provides 
the required distribution of the flow variables, \p, w, H, etc. 
in the physical space. The velocity components W^, in the 
physical space can be related to the transformed ^ derivatives 
using equations (9) and (31) with the following result: 


W = 
m 


M , 3x 3i{> 3x 3jK , 
bpr '35 3n ~ 3n 35 J/ 


w _-_m_ f ii li _ li ii) /t 

“ bpr l 3n 35 35 3n J/ 


(37) 
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Boundary Conditions in the Transformed Domain 

The boundary conditions for equations (32) to (34) are 
obtained by a proper transformation of the boundary conditions 
stated earlier in Section 1 to the new system of coordinate (5,n) • 
This procedure yields the following relations. 


a. The Upstream Boundary AN 


M = _i fll _ M tanR) 

3 5 2ir/Z l 35 3 5 inlet 


H _ 

an 


w “ - (2Si sina) inlet 


H ^ C p T total ^ inlet 


(38a) 


(38b) 

(38c) 


b. The Periodic Flow Boundaries AB, NM and FG, IH : 

Figure 6 shows the grid ordering system in the transformed 
domain. The grid rows along j = 1 and j = p correspond to the 
circumferential boundary of the blade-to-blade passage AB, NM, 
and FG, IH. The grid rows along j = - 1 and j = p + 1 are 
exterior to the blade rows and are reserved in the computational 
procedure for the enforcement of the periodicity condition. 
Accordingly, along j = 1 and j = p, the following conditions apply 


+ ±,-1 = ^p-l" 1 ' *i,l = *i,p 
“i,-l = w i,p-l ' “i,l = “i,p 
H i,-1 ~ H i,p-1 ' H i,l ~ H i,p 


1 , 

*i,2 

^i,p+l ^ 

(39a) 

r 

“i,2 

= “i,P+l 

(39b) 

i 

H i,2 

= H. 

1,P+1 

(39c) 
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where i, is the finite difference grid counter in the 5 direction. 
Similar relations are used for the grid points along FG and XH. 

c . The Blade Surfaces MI and BF : 

Laminar flow case 


iff = 0 

(Along MI) 

(40a) 

♦ = 1 

(Along BF) 

(40b) 

% - ■ 

[a( Vi'V + ! “w+i +Bl ' c 

(40c) 


(For both Mi, BF) 



H = C t + “ r (for both MI, BF) (40d) 

w p w £ 
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M 
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In the a b ove equations the subscript (w) denotes a blade boundary 
point and the subscript (w+1) denotes a point in the flow field at 
a distance An away from w. The vorticity boundary condition given 
by equation (40c) is derived using a Taylor's series expansion 
for the stream function about the blade surface, the boundary 
conditions given by equation (23) and the stream function 
equation (33) . The formulation of the vorticity boundary condition 
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as given by this equation is second order accurate and allows 
implicit treatment in the solution technique. 

d. The Downstream Boundary GHi 


8# 1 

( H - If fc^axit 

(41a) 

3£ 2j/Z 

to ~ — [252 

sino) e X it 

(41b) 

|l-o 

3? 

' 
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iv . Pressure Distribution 

In the stream function-vorticity formulation of the flow 
equation, the pressure does not appear explicitly in the problem. 
Therefore, indirect methods must oe used to evaluate the pressure 
distribution through the flow field. An accurate approach 
consists of forming a Poisson equation for the pressure using 
equations (7) and (8) . The Poisson equation is then solved to 
determine the pressure, subject ~o Neumann -type boundary conditions 
provided by the momentum equations. A variety of numerical 
procedures for this type of solution are given in reference [13] . 

The principal problem encountered with such an approach is that 
the Neumann boundary conditions should be formulated carefully 
such that they are compatible with the source term in Poisson 
equation. Because of the trun<_ ~tion errors, the boundary values 
fail usually to meet this constraint, resulting in a slow divergence 
of the numerical solution. An alternative approach to obtain a 
pressure distribution is to perform a line integral for the 
pressure gradients |£, -|E along a contour in the flow field. 

This approach usually yields reasonably accurate values for the 
pressure along the surface of a smooth body without sharp corners 
as is the case in the present flow configuration. 

The pressure gradients td -|^ are related to the velocity 
components W^, W^, the vorticity j and their derivatives by the 
following relations: 


28 


and 
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vi i 3A 3W 
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34> 3W m 3x 3W 4. 
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In deriving these relations/ use has been made of the momentum 
equations (7) and (8) , together with the transformation rela- 
tions (lb) through (10b ) , as well as the definition of u 
given by equation (14) . . 

Once a solution to the transformed flow equations has been 
obtained/ the right hand side of equation (42a) and (42b) can 
be evaluated using second order central difference for the 
5 and n derivatives. Prom a mathematical view point/ equations 
(42a) and (42b) represent algebraic expressions for the various 
pressure gradients. The pressure distribution over the whole 
field can therefore be obtained by performing a line integral 
for these algebraic expressions. Starting from the upstream 
inlet boundary (AN) of Figure 6 where the pressure level is 
known, the pressure distribution along the mid channel line (L-L) 
is obtained using the following integration formula: 
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(43) 


h I! 35 =-Pi-i + I 


A similar expression is used to integrate equation (43) in the n 
direction. In this case the integration starts with the predicted 
pressure values along the mid channel line (L-L) and proceeds 
towards both the suction and pressure surface of the blades. 


3 . SOLUTION PROCEDURE 


The flow governing equations introduced earlier as equations 
(32) through (34) comprise a system of coupled nonlinear elliptic 
partial differential equations that must be solved, subject to 
the boundary conditions (38) - (41) , to provide the details of 
the flow pattern within the blade-to-blade channels. Since the 
flow equations are not tractable to analytical solutions, a 
numerical solution which is based on the finite difference method 
is used. The following simplifications are made in order to 
reduce the complexity of handling the numerics of the problem. 

First, the flow is considered to be incompressible. This decouples 
the energy equation from the momentum equations. Thus, it is 
possible to solve the i|» and w equations to obtain the velocity 
distribution which are then used in the energy equation to determine 
the enthalpy distribution throughout the flow field. This 
assumption is mandatory in order to develop a method of solution 
which could be later expanded to take into account the density 
variation within the flow. Second, a non -contracting body fitted 
coordinate system is employed during the numerical solutions. 
Consequently, all first order terms containing the transformation 
parameters a, x,in equations (32) - (34) vanished as pointed out 
in Section 2 . It should also be remarked that using non- 
contracting body fitted coordinates implies that the spacing 
of the coordinate lines, in A direction, in the physical domain, 
is uniform. 
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Before writing the governing equations in finite difference 
form, it is found more convenient to express the flow variables in 
the following dimensionless form: 


* 

to = 


air, 


( V 0 


Re 


* 

H 


H 


"ft, 









(44) 


Where the subscripts o, t denote the condition at the upstream 
boundaries (A-N) and M-B) respectively (see Pig. 3) . Using the 
above equation, one can express the flow governing equations in 
nondimensional for as follows: 


Stream Function Equation 
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Vorticity Transport Equation 
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Energy Equation 
* 


b* 
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; ** p 3n ^Re Pr an ; ‘ Y an l Re Pr 3n J 


(47) 


where z is the number of blades, D and G 4 are the nondimensional 
equivalent of the source terms D, G, in equations (34a) and (34b) . 


Numerical Solution 

The derivations of the finite difference equations is followed 
by a description of the iterative procedure used to obtain the 
numerical solution. 

The Finite Difference Equations 

A noteworthy comment is to be made before expressing equations 
(45) - (47) totally in a finite difference form. This relates to 
the effective handling of the first order convective terms in 
the above equations. As pointed out by Roache [13] , these terms 
can destroy the diagonal dominance of the matrix of the finite 
difference equations to be solved at high Reynolds numbers. This 
in turn causes inversion instabilities that produce spatial 
oscillation "wiggles” in the final solution for the flow variables. 
To eliminate these instabilities, a windward difference technique 
is used to model the longitudinal convective terms. Thus, when 
the local value of the £ component of velocity is positive, the 
convective terms including or -Iv— are evaluated with a 

backward difference. 


ae ae 

On the other hand, when the velocity is 


negative, a forward difference approximation is used for either 
or -If^- . In a similar fashion, the windward differencing 
technique is used to control the difference representation of the 


normal convective terms containing and 

(46) and (47). 


3H* 

3n 


in equations 
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Referring now to Pig. 6, it can be shown that the 
finite difference representation of equations (15) - (47) may 
be written as 


♦±,j = & 1 + A 2 + A 3 ♦i.j+l + a 4 *i,j-X 


+ A r 


( 48 ) 


and 


■ft ic it ^ ^ 

tu. . = B-. tu . . + B_ (I). _ t* B- to. . - + B. tu. . 1 

1 f ^ JL 1 "T" 1 ^ £» X “ i / _ 3 1 ^ J 4* 1 4 ^ ] "*X 


+ Bh 


(49) 


and 


it •k * it it 

H. , » C- H. . , . + Ch H. . • + C_ H. . t . + C. H. . , 
i/j i 1+1,3 2 i-l,3 3 l/3+l 4 i/3-l 


+ C r 


(50) 


Where the coefficients A^/ B,. and C^, in the above equations 
serve as the explicitly known source terms. The value of these 
coefficients as well as the cofficients A^ through A a , B^ 
through , and D n through / are given in Table 2 . 

Iterative Procedure 

For the numerical solution of our system of coupled finite 
difference equations (48) - (50) , a wide variety of classical 
matrix iteration techniques are available. The point SOR method 
of solution has generally been acknowledged as being the best 
of these iterative techniques/ because of its effectiveness 
and simplicity of application. Within the category of the point 
iterative method/ there is a choice of alternate methods. 

In the present report/ the Gauss-Seidel method is used, for it 
is known to yield more rapid convergence solution and places less 
demands on computer storage. 
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TABLE 2, COEFFICIENTS OF THE FINITE DIFFERENCE EQUATIONS 


A i = 26 An/A?/A o /(b * +1 + b *) 
A 2 = 2SAn/A£/A o /(b*_ 1 + b*) 


A 3 = y AC/An/A Q /b i 


U) 


A 4 = Y A^/An/A Q /b i 


*"t ^4- * 2 2 

Ag = A? An [-^ ^ w PTCT z/2ir 
o o 


a ^ /i ^ 1 k i •> q 3 /I 3 ^ > -1 / n 

£ W ^ 9n } " £ 3n 3/ o 


where 


A — A n + A„ + A-. + A. 
0 12 3 4 


t All unsubsuript quantities in this table are evaluated at the grid 
node (i,j). Note that the stream channel thickness b* is constant 
for all j . 1 
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°21 " J b o ^*22 + 1*22 1 > /8b * 

CL, = 6 An/ A£/Pr/ [ (Re) . . + (Re). .] 
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* The Prandtl Number is assumed to be constant in this table 


In using the Gauss-Seidel method, two cycles of iterations 
are required to obtain a complete solution to equations (48) - (50) . 
The first cycle consists of the iterative solution of the stream 
function and vorticity equations . When a converged solution is 
reached for these two equations, the energy equation is then 
solved in the final cycle. During the first cycle of iterations, 
a successive over-relaxation, using an optimum over-relaxation 
factor, K.£, is used. The iterative procedure is given by 

0 = R f e n+1 + (1 - R f ) 8 n (51) 

where 9 denotes a general flow property and includes ^ and to . 

The superscripts (n+1) , and n represent an iteration counter. 

For the final cycle of iterations it is found that using 
an under relaxation factor, U, improved the convergence character- 
istics. The under relaxation relation used is 

H = 0 H n+1 + (1-U) H n (52) 

In the present solution the optimum values for the different 
relaxation factors are determined by trial and error. 

Boundary Conditions in the Iterative Scheme 

The general recursive formulae for the iterative solution 
as given by equations (48) - (50) are only applied at the interior 
nodes of the flow domain. Near the various boundaries equations 
similar to (38) - (41) are used. Along the upstream and downstream 
boundaries, equations (38) and (41) are introduced in the solution 
procedure in a straightforward manner with minimum compu- 
tational effort. However, along the periodic boundaries AB, NM, 

FG and IH, the boundary conditions provided by equation (39) 
have a unique feature. Specifically, neither the boundary values 
nor the derivatives of any of the flow properties are specified 
apriori along these boundaries. The only information supplied 
by equation (39) are the equality of the function values and 
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their derivatives at every two corresponding points. Therefore, 
a special approach is required to determine the values of the 
dependent variables if /, to, and H, along the periodic boundaries. 

The approach to be taken here relies on a modified evaluation 
of the governing equation for each flow property to determine 
the boundary condition for the flow property in question. The 
modification is introduced in an attempt to insure that equation 
(39) holds, that is to satisfy the periodicity condition. The 
numerical procedure is as follows. 

Let the generalized finite-difference form of any of the flow 
governing equations (45) - (47) be represented by, (see Fig. 6 
for subscript notation) : 




= a 


1 


+ a„ 0 


2 i-1, j 


+ a. 


0 . . , , + 

i/3+l 


a 4 


+ a 5 (53) 

Where 0 denotes a general flow property (ill, to, and H) . The 
coefficients a^ through a^ depend on the particular equation 
used, and are obtained by a proper combination of the coeffi- 
cients A^ through Cg given in Table 2. For the grid points 
along the boundary AB of Fig. 6, the above equation is modified, 
using the periodic boundary conditions given by equation (39) 
to the following form: 


0. . = a. 0 . . _ . + a„ 0 . t . +a„ 0. „ + a. 9. . 

1,1 1 1+1,1 2 1-1,1 3 1,2 4 i,p-l 


+ a 5 + X 1 


(54) 


where x-j_ = 0 for to and H; 

x x - - a 4 f °r 4> • 


(54a) 
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i a, ■■ i 


The points along NM are not part of the solution regions, since 
the value of the dependent variables at each of them is just 
equal to the corresponding point along AB. The equation for 
the first mesh line below NM must be modified by substituting 
the periodicity condition given by equation (39) into equation (53) . 


0 i,p-l “ a l 4 * * * * 9 i+l,p-l + a 2 0 i-l,p-l + a 3 0 i,l + a 4 0 i,p-2 


+ a 5 + X 2 


(55) 


where x 2 = 0 for “ ant ^ H; 

X 2 = a 3 for i|> . (55a) 

A similar approach can be applied along the other boundaries 
(FG and IH) where a periodic condition exists. 


4. RESULTS AND DISCUSSION 


The equations formulated in Section 2 are programmed for 
numerical solution using the finite difference technique discussed 
in Section 3. The program is arranged to handle general flow 
within turbomachinery, which may be of the axial, radial or 
mixed flow type. In general, the program requires as an input 
the configuration of the stream channel annulus S^, the inlet 

flow conditions, the rotational speed of the machine, and the 
blade geometry. Recalling that all the flow calculations are 
cairried out in the unit square of the transformed domain, there- 

fore the blade input geometry is supplied to the program in the 

form of the transformation parameters 6, 0, y and J. As pointed 
out in Section 2, these parameters may be specified for any blade 

geometry using Thompson code for the automatic numerical generation 
of boundary fitted coordinate system [12 ] . 


40 


The program output consists of the distribution of the 
stream function, the vorticity, and the static pressures within 
the blade passages. The variation of meridional and tangential 
velocity components from blade-to-blade and from the inlet of 
the machine to its exit are also generated. In cases where 
blade cooling is considered, the program has the capability to 
generate the temperature distribution within the blade passages. 

In order to keep the computer time within reasonable limits, 
(usually less than 5 minutes on an AMDAHL 470) , the flow domain 
has been divided, for all calculations, into 30 step sizes in n 
direction and 40 in the i; direction, with the greater number of 
nodes distributed in the meridional direction. 

Five flow cases are investigated using the developed program. 
The main purpose was to check the accuracy of the present method 
of analysis in predicting the actual flow behavior within 
turbomachine channels . The accuracy of the method was confirmed 
by a comparison with available experimental data. Four of those 
investigated cases were concerned with inward flow situations, 
while the fifth one dealt with an outward flow case. In all 
cases investigated, the flow was considered to be incompressible 
and having a constant effective viscosity, u e * The blade surfaces 
are assumed to be adiabatic with no heat sources or sinks. 

Inward Flow Cases 

These flow cases are those of a radial inflow turbine whose 
rotor consists of eight radial straight blades. A full description 
of the rotor geometry is given in Fig. 7. The primary reason for 
the selection of this specific rotor is that a substantial 
amount of experimental data is available for It (reference [20]). 
Thus, besides providing a basis for comparison with the theoretical 
predictions, the experimental evidence is used to show avenues for 
future development in the present method of analysis. 

The flow patterns are investigated on a blade-to-blade stream 
channel, S^, located midway the passage depth of the rotor as 
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shown, in Fig. 7. The results are presented over a wide range 
of operating conditions, which are summarized in Table 3. The 
coordinate system used in the solution is illustrated in Fig. 8. 
The results presented include stream function contour plots, 
velocity profiles across the rotor passages together with some 
information concerning the pressure distribution within these 
passages. 

The streamline contours for the four inward flow cases of 
Table 3 are shown in Figs. 9a and 9b. The streamlines are plotted 
for the region between a pair of blades, represented by the 
heavy thick lines . The streamlines are designated by a stream 
function ratio ty/ty to tal suc ^ the value on a streamline 

indicates the ratio of flow through the passage between the 
streamline and the pressure surface of the blade. Thus, for the 
given channel configuration, the streamline spacing is indicative 
of the velocity relative to the rotor, with close spacing indi- 
cating high velocities and wide spacing indicating low velocities. 
For the operating conditions corresponding to case 1, as shown 
in Fig. 9a, it is observed that a recirculating 'ddy began to 
form near the pressure surface of the blades. As the rotating 
speed increases, the recirculating zone grows much larger, as 
shown for case 3 in Fig. 9b. The relative velocity near the 
suction side of the blades increases in the later case. This 
effect may be attributed to the fact that the effective sectional 
area of the rotor decreases with the growth of the recirculating 
zone. Since large recirculatory zones cause higher losses in 
total pressure, it is desirable to avoid them through efficient 
rotor design and proper selection of the oeprating conditions. 

From the inspection of Fig. 9, it can be concluded that the size 
of the recirculating eddy depends upon the relative magnitude of 
the flow rate (M/T^/p^) and on the rotor speed, N//T^. These 
zones generally can be reduced by increasing the mass flow rate 
and/or decreasing the rotor tip speed. 

The most remarkable feature of the present results is the 
good agreement obtained between the predicted flow behavior and 
the experimental evidence taken from reference [20]. In all 
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cases studied, the size and the extent to which the recirculating 
zone grows compares favorably well with the experimental data, 
as indicated in Fig. 9. 

Figures 10a and, 10b show the predicted meridional velocity 
distribution across the blade passages at three different radial 
locations. These locations are selected to correspond to a 
radius of 23.2, 15 and 6.8 eras, respectively. The dotted line 
shown on each figure represents the velocity distribution for 
an inviscid solution. The results of the viscous flow analysis 
show a large variation in the meridional velocity profiles as 
the flow travels downstream. The profile distribution, at 
stations located away from the turbine inlet, indicate that 
regions of high meridional velocities are shifted towards the 
blades suction surface as shown in Figs. 10a and 10b. While, 
regions of high relative meridional velocities are observed to 
exist near the blades pressure surface at subsequent downstream 
stations. Compared with the viscous calculations, the inviscid 
flow solution predicts a completely different flow behavior, in 
this respect. Moreover, in some flow situations where severe 
changes take place near the rotor tip as in case 4 in Fig. 10b, 
the inviscid flow solution fails completely even in predicting 
the flow characteristics. All these factors, in addition to the 
existance of reversed flow regions near the blades pressure 
surface make it clear why the inviscid solutions always fail to 
produce a realistic prediction of boundary layer characteristic 
parameters in rotating machines when used in conjunction with 
standard boundary layer analysis. 

Figures 11a, b,c show the pressure variation between blades at 
four radial locations corresponding to radius, r = 13.1, 16.3, 

19.5 and 22.7 cms. The static pressure, p, is plotted in these 


figures using the nondimensional quantity. 



(P x -P) / 


p 1& r 


2„2 
tin 


2g 


where p^, are the mean static pressure and density at rotor 
inlet respectively. The experimental measurements of reference [20] 


43 


are also shown in the same figures. As with the stream function 
and velocity results, the present method of analysis provides a 
good prediction of the pressure distribution over a wide range 
of operating conditions. On the whole, the value of C is 
observed to be larger on the suction . side of the blades and 
decreases near the pressure surface, with smaller changes in C 
values at higher rotational speeds. In the region lying between 
the center of 'the passage and the pressure surface, the 
values become smaller as the rotational speeds decrease. These 
observations imply that the static pressure drop from the rotor 
inlet increases with increased rotational speed. Near the 
suction surface, the pressure drop increases with the decrease 
of rotational speed. Such a tendency is remarkable, particularly 
near the rotor inlet. 

Some comments might also be made concerning the discrepancies 
observed between the predicted and experimental values of pressure 
distribution in case 4 near the rotor exit (i.e. at r = 13.1) . As 
reported in the experimental work of reference [20] , the operating 

A 

conditions for this case cause the flow to be heavily separated . 
The existance of large separation zone within the rotor is 
believed to modify the channel shape in such a way as if the 
zone acts as a pseudo blade. This in turn affects the pressure 
distribution in the manner shown in Pig. 11c. Therefore, it 
appears that in order to deal with flow cases where heavy 
separation is encountered within the rotor channel, one has to 
incorporate a zonal model for such separation in the present method 
of analysis. Such development in the method of solution is 
undoubtably essential for further use of this method of analysis 
in aerodynamic improvement and performance prediction of 
turbomachines . 


t It is to be noted that this type of separation is not a regular 
two dimensional one, but rather a three dimensional type. 
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Flow in Radial Compressor 

The capability of the present method to analyze the flow in 
diffusion cascades is examined by studying the flow behavior 
within the rotor of a radial bladed compressor. The rotor 
profile is shown in Fig. 12. The rotor has 12 straight radial 
blades . The blade-to-blade shape in the physical domain as well 
as the coordinate system used in the solution are shown in 
Fig. 13. Additional summary data for the solution appear in 
Table 3 (as case number 5) . A typical distribution of the flow 
properties on the blade-to-blade computational surface, S^, of 
Fig. 12 are calculated and the corresponding results are cal- 
culated and the corresponding results are presented in Figs. 14 
and IS. 

Figure 14 shows a comparison between the predicted stream- 
line contours and those determined experimentally in reference [21] . 
The experimental evidence was obtained by tracing photographs of 
streak lines from the rotor segments under the same operating 
conditions reported here. Good agreement is generally observed 
and it should be particularly noted that the shape and the size 
of the recirculating eddy compare favorably well in both cases. 

The predicted meridional velocity profiles across the rotor 
passage at three different radial locations are illustrated in 
Fig. 15. Shown also in the same figure, the calculated velocity 
distribution using an inviscid blade-to-blade analysis. The 
inviscid solution, although showing the existance of negative flow 
regions as exemplified in Fig. 15b, over estimates the size 
of the recirculating eddy noted in Fig. 14. This overestimation 
is supported by the existance of large negative meridional 
velocities near the blade suction surface. In an actual case, 
boundary layer phenomena are expected to reduce the effective 
flow area of the passage, thus increasing the volume flow rate 
per unit area through the effective area and thereby reducing 
the size of the eddy. This is exactly the same result obtained 
using the present viscous flow solution. 
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TABLE 3. PARAMETER FOR THE NUMERICAL SOLUTIONS' 
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5 
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In all cases studied 

radius ratio r/r 

radius ratio r/r^. P 

tio 


, the upstream boundary AN of Figures 8 and 13 is located at 
of 1.35, while the downstream boundary HG is located at 
of 0.254. 


** H designates the rotational speed (r.p.m.) 


t For all cases, the blade surfaces are considered to be adiabatic and the flow ls 
assumed to be laminar with constant viscosity, y g . 
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On the whole, the present method of analysis provides a 
good prediction of the actual flow behavior within the passage 
of turbomachine rotors. The preservation of the ellipticity 
of the flow problem is believed to be the major element that 
results in such good prediction. The ellipticity is preserved 
through the consideration of all the diffusion terms of the 
governing equations during the solution procedures. The method 
of analysis proves also to be of acceptable accuracy and 
provides invaluable information on the rotor flow characteristics. 
This is evidenced by the good agreement obtained between the 
predicted results and the available experimental data. 
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NOMENCLATURE 


A 1' A 2' A 5 

B 1 ,B 2' B 5 


C 1' C 2' C 5 

G 1' G 2' C D 


D 


m 


(e n>? 

( ®n*n 

(i t ) c 

(i t>n 

E 


G 1' G 2' G 3' 

G 4' G 5 


/ h^ r 
H 
J 

l 

\ 

m 

M 


n 


Coefficients of the finite difference equations. 
Coefficients of the finite difference equations . 

Normal stream annulus thickness , m. 

Coefficients of the finite difference equations. 

Constants in the turbulence model. 

Specific heat at constant pressure, J/ (Kg) (°K) . 

Dissipation function. 

Unit vector in m-direction. 

Unit vector in <j>-direction. 

Unit vector normal to a constant line. 

Unit vector normal to a constant n-line. 

Unit vector tangent to a constant £-line. 

Unit vector tangent to* a constant rt-line. 

Kinetic energy of turbulence, J/Kg. 

Donoting source terms in the flow governing equations. 
Static enthalpy, J/Kg. 

Scale factor for the orthogonal curvilinear "oordinates . 
Total enthalpy, J/Kg. 

Jacobian matrix, E. (B.2). 

Mixing length, m. 

Meridional distance, m. 

Mass flow per blade flowing through the stream 
annulus. Kg/sec. 

Outward unit normal to the stream surface, 8^, see Fig. 2. 
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N 


N 1' N 2 

AN 


P 

Pr 

Pr 

r 

R 


Re 



t 

T 

V 

V 

w 

w 


W W 3 


WL 


m 


w. 


x 


W X 3 


Outward unit normal to blade surface. 

Components of N in m, <f> directions respectively. 

Distance of the near-wall grid point (w+1) from the 
blade surface , see Pig. 5. 

2 

Static pressure, N/m or blade pitch. 

Effective turbulent Prandtl number. 

Molecular Prandtl number. 

Radius from axis of rotation, m. 

Universal gas constant, U/(Kg)(°K). 

Reynolds number. 

Schmidt number for kinetic energy of turbulence. 

Schmidt number for dissipation of kinetic energy of 
turbulence. 

Time. 

Temperature, °K. 

Absolute velocity vector, m/sec. 

Magnitude of V, m/sec. 

Relative velocity vector, m/sec. 

Magnitude of W, m/sec. 

Components of W in x^, and x^ directions, 
respectively. 

Meridional component of the relative velocity vector, 
m/sec . 

Tangential component of the relative velocity vector, 
m/sec. 

Stretched meridional coordinate, Eg. (26) . 

General orthogonal curvilinear coordinate. 
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z 

Z 

a 

3 


Axial coordinate, m. 

Number of blades. 

Angle between m and z, rad., see Fig. 2. 

Coordinate transformation parameter, Eq. (29) , 
or angle between relative velocity vector and 
meridional plane, rad., see Fig. 2. 


Y 

Coordinate 

tr an s formation 

parameter , 

Eq. 

(29) . 

6 

Coordinate 

transformation 

parameter, 

Eq. 

(29) . 

P 

3 

.Fluid density, Kg/m . 




T 

Coordinate 

transformation 

parameter. 

Eq. 

(32b) . 

a 

Coordinate 

transformation 

parameter , 

Eq. 

(32b) . 


V 

u e 

S2 

0 

<P 

* 

CO 

e 

n 

At 

Superscripts 


n 


Kinematic viscosity, m /sec. 

2 

Effective viscosity, m /sec. 

. 2 

Laminar viscosity, m /sec. 

Rotational speed, rad/sec. 

Denotes general flow property and includes iff, co, and H. 
Relative angular coordinate, rad., see Fig. 2. 

Stream function. 

Vorticity, 1/sec. 

Dissipation of Kinetic energy of turbulence, J/Kg. 
Boundary fitted coordinate, see Fig. S. 

Boundary fitted coordinate, see Fig. 6 . 

Time step counter. 

Mean value 

Denotes nondimensional quantity. 

Iteration counter. 
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Subscripts 

e 

l 

m 

w 

w+1 

* 

if j 

inlet 

exit 

total 

tip 


Effective 
Laminar . 

Meridional component. 

\ 

Wall value. 

Pertaining to points in the flow field at a 
distance An away from w, see Pig. 6. 

Tangential component. 

Denotes field poistion in (£,n) domain, see F.ig. 
Inlet or upstream. 

Exit or downstream. 

Total conditions. 

Rotor tip. 
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STREAM SURFACES CONTOURS, OBTAINED FROM A 
MERIDIONAL FLOW ANALYSIS, FOR A MIXED FLOW 

E. 



Grid Pdinla * 




Detail A-At Grid Pointo Hear the Ulade 
Surface , 



Detail B-Bi The Hear Hall Region* 


FIG, 4. GRID POINTS NEAR SOLID BOUNDARIES, 












Angular Coordinate, <p radians 











MERIDIONAL VELOCITY DISTRIBUTION, Tf_/U 


VISCOUS 


INVISCID 


Angular Coordinate, 
(f> , degrees 




PREDICTION CASE 1 


Angular Coordinate 
<j>, degrees 
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PREDICTION CASE 2 


FIG, 10a, NONDIMENSIONAL VELOCITY DISTRIBUTION AT 
DIFFERENT RADIAL LOCATIONS, 


VISCOUS 


Angular Coordinate, 
' <f>, degrees 



r — 23 . 2 cm 
= 39.3 m/s 


1 0.4 



INVISCID 

Angular Coordinate, 
4> , degrees. 



PREDICTION Ch-'.l 3 


PREDICTION CASE 4 


FIG, 10b. NONDIMENS IONAL VELOCITY DISTRIBUTION AT 
DIFFERENT RADIAL LOCATIONS. 
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Angular Coordinate, (ji, degrees 

EXPERIMENTAL CASE 2 


FIG. 11a. COMPARISON OF PREDICTED STATIC PRESSURE DISTRIBUTION WITH 
EXPERIMENTAL DATA OF REFERENCE (21) . 
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CASE 5 


COMPARISON OF PREDICTED STATIC PRESSURE 
DISTRIBUTION WITH EXPERIMENTAL DATA OF 
REFERENCE (21), 
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a) HERIUIGHM, VIEW b) niADE-TO-ULA.DK VIEW 


FIG. 12. HUB-SHROUD PROFILE WITH THE STREAM SURFACE, Sj, USED FOR 
BLADE-TO-BLADE ANALYSIS (COMPRESSOR CASE). 


Angular Coordinate, <}> radians 











MERIDIONAL VELOCITY DISTRIBUTION, W/U 



Angular Coordinate, <p, 


FIG, 10. VELOCITY PROFILES ACROSS THE 
PASSAGES AT DIFFERENT RADIAL 
(CASE 5 IN TABLE I), 




APPENDIX A 


DERIVATION OF STREAM SURFACE EQUATIONS 


The detailed derivation of the equations governing the 
fluid flow on the blade-to-blade computational surface r S^, 
of Eig. 2 is presented in this appendix. The complete equations 
of motion that describe the flow behavior in a turbomachine 
passage, are first written for a general orthogonal curvilinear 
coordinate system. A transformation of the resulting equations 
'to the computational coordinates system is then outlined. 

Principle Equations in Orthogonal Curvilinear Coordinate 

For a general orthogonal curvilinear coordinate system 
with scale factors h», h 2 , h 3 , the Navier Stokes 
equations given as equations (1) and (2) in Section 1 may be 
written in the following form: 

Conservation of Mass: 


3p . 1 3 . , 1 3tp 'V 

at + iq asq lpW x ) + 5^ ax 2 


* hr 2 < pw 3> 


PW ] 

pW. 


* h. 



°“2 

. 1 


5x l h 3 

ax i 


3h l . i_ 

3h 2 

3x 3 

h 2 

3x 3 ' 


pw_ , 3h 

* + hT t 35 


) = 0 


2 h 3 3x 2 


(A. 1) 


where W are the components of the velocity vector W 
in the x^, X£ and x 3 directions respectively. 
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Conservation of Momentum in xt Direction 


P [ 


^2 ”2 3h 2 

' h 2 3x 2 h 3 3x 3 _ h l h 2 3x l 


at 


h^ 3x^ 


W. 3h ’ 

— — — 1 

h 2 ax 2 J 


W 3 , W 1 ^^1 W 3 ^^3 o 2 1 a r 2 

(r - *=* - 57 asf> + 2 («2 W 7 - <W" 3- C- IS-] 


“1 3 

1 9 P 

3x^ 


h„ 9x. 


‘2"3 


'3 2 


h i 3x i 




ll l^ 1 2 il 3 


3x n 


1 3W 1 

[2u e h 2 h 3 + 


h l 3X 1 h l h 2 3x 2 


w„ ah. w ah. 

2 1 + r~^~ TTri) 1 


h 3 h 1 ax 3 


3x, 


[ « e h i h 3 <hf k 


w. 


a 


w. 




a 

3 x 3 

tu e h 

l h 2 

h 3 3x 3 

£ 

h l 


\ 

a 

W 2 

n 2 

h l 

h l h 2 

3x i 

h 2 

»e 

h i 

{-A 

h 3 

3 

W. 

<4> * 

h. \ 

h^h 3 

ax 3 

ll ^ < 

a) 

<N 

(i- 

h 2 

aw 2 


3h 2 

h l h 2 

3 x 2 

+ 

h 2 h 3 

3X 3 

**• 

t±- 

h 3 

aw 3 


3h 3 

h 1 h 3 

3 x 3 

h 3 h l 

3X 1 


3 a 


w. 


h i 3x i h 3 


2_ (” 1 , } . !^1 EODDCIBILITY of the 

3x 2 V ' 


— — — nvi/ 

3x 2 JGINAL page IS POOP 


w ah, 
( — =L) } • ±. 

ax, h 3 ; 3 x 3 


W 1 3h 2, 


h X h 2 3X 1 


w. 


3h. 


h 2 h 3 3x 2 ' 


3h. 

t 

9x, 


ah, 

aiT! 


(A. 2) 


where 

7 *W = 


1 

" h l 

aw x 

3x i 

. 1 
" h 2 

aw 2 

ax" 

*k 

aw 3 

3x 3 + 

W I 

h l 

(i. 

h 2 

^2 

3x i 

+ j“ 

h 3 

ah. 

ax i 

W 2 
+ h" 
h 2 

(L. 

h i 

3h l. . 

3x 2 

*k 

ah, 

3x 2> 

W 3 
+ IT 1 
h 3 

<jj- 

h l 

3h l 

ax 3 

h 2 

3h, 
3x ' 

J 

(A. 2a) 
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Conservation of Momentum in Direction: 


aw. 

* i 

: at 


o + Ha * Hi Oh * h. ^ ”-1 ,h. Sh 3 w a 3h 2, 
h l 3x l h 2 3x 2 t h 3 3x 3 h 2 ( h 3 ax 2 - h 3 3x 3 > 


£. f^A ^2 _ j_ ^, n „ n r - \ n 2 3r 2 7 

' Ct - h, 3xJ + 2(J2 3 W 1 “ *W~ 2hJ 


h n h., 3x. 


■ EJ llj ■ f EJ lw e 7-N) 


w_ 


w. 


h l h 2 h 3 3x l 


f a CTi h .. r_2 3 ( _ 2 s . ^ 3 "1.., 

"•av [n«h.,jj (. a »v, ) 1 , a,, ft, ))1 


2 3 e l h x 3x x 'h 2 


h 2 3 x 2 


+ 3 

[2U e h l h 3 

(i_ 

h 2 

3W 2 

+ W 3 Sh 2 + 

W x 3h 

3x 2 

3X 2 

h 2 h 3 3x 3 

h l h 2 3X 

, 3 

--5S7 

[l, e h l h 2 

,5a 

l h. 

a 

3x. 

(% + ^ 3 
'h. J h_ 3x_ 

W 2 

(“))]} 
XI « 


)] 


+ r 


t^ 3 


W 3 h 2 3 

Cc 1 ) + 23 


W 2 . 3h 2 

h 2 h 3 l, h 2 3 x 2 'h 3 v ‘ h 3 3 x 3 ^ 2 ^ 3x 3 

M e f ^2 3 *2 , 3 ^1 . ^2 

h 2 h l h l 3x l h 2 h 2 3x 2 h l 3x l 


2 “e 

% 

aw 3 

W 1 

-f- 

Sh 3 

+ W 2 

3h 3,. 

Sh 3 

li 2 h 3 

3x 3 

h 3 h l 

3x l 

h 2 h 3 

3X 2 

3x 2 

2V e 

l k 

3w i 

t W 2 

3h l 

+ W 3 

3h, 

1 

3h l 

h 2 h l 

3x X 

h l h 2 

3x 2 

h 3 h l 

3x 3 

3x 2 


(A. 3) 
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RE ?RODUClBlL ITr 

Conservation of Momentum in x 3 Direction ; ~ 


n-\ ... 


3W, W. 3W, W, 3W, W, 3W, W, W, 3h_ W, 3h, 

0 r. ± -r — - 4. -A 3 + __3 3 1 /_1 1 3 J} 

ML 3t Sx^^ h 2 3 x 2 h 3 3 x 3 h 3 l h 1 3x 3 N h^ 3x^ 

w, w 1 ah3_w i 3h i di_ar* 

h 3 ( h 2 9 x 2 h 2 ax 3 ’ T 2tB l W 2 n 2 W l )_ 2 h 3 SXj 1 


- — — S — - — — (u v -W) 

h 3 3X 3 3 h 3 3 x 3 ' y e V WJ 


ft, a W. ft- » W- 

h - Ch-h-U A f— (r-i) + ^ f— (r 2 ))] 

h l h 2 h, 3 3x l 2 o e h 3 3 x 3 h-j^ h i 3x x h 3 


V, 


a ft, - W, ft, a W, 

* 3xJ [h l h 3 y e ( hj 3x^ ( ftj } + hj 3^ ( hJ J } 1 


3x, 


[2h.,h n U (: 




W, 


h l h 3 Lh 3 3x 3 V 


3W 3 

W 1 

3h, 

W 2 

3x 3 "** 

h 3 h l 

3x i ’ 

h 2 ft 3 

ft- 


W- 

3h_ 

+ *7 

3 

3x^ 

( ir n * 

3 

3x., 


L )]> 


h 2 h 3 


ft, , W h - w , 3ft a 

r_3 3 /_!) + _1 3 (_2) 1 . 3 

h 2 3 x 2 h 3 ft 3 3 x 3 h 2 3 x 2 


2u 


h l h 3 


x 3W 1 w 2 w 3 3h 1 

h l 3x l h l h 2 3x 2 h 3 il l oX 3 


1 


3h ] 

3x, 


!^s_ f i_ !!a , J!i_ ^2 + J0i_ ^ ^ 

h 2 h 3 h 2 3x 2 ' h 2 h 3 3x 3 h l h 2 3x l 3x 3 


(A. 4) 


In the above equations £2^, Q 2 , £2 3 are the components of the rotor 
angular velocity, Q, in the x^, x 2 / x 3 directions respectively. 


76 


Flow Equations on Blade-to-Blade Computational Surface Si: 


Starting from equations (A.l) , (A. 2) and {A. 3) , the flow 
equations over a blade-to-blade computational surface is formu- 
lated. In deriving the required equations, one has to describe 
the computational surface in terms of the orthogonal curvi- 
linear coordinates x^, x 2 and x^, 

S 1 (x 1 , x 2 , x 3 ) = 0 (A. 5) 


Equation (A. 5) is used to relate any flow property q of the three 
dimensional flow field with the same flow property q on the 
surface S^. In general 

q a x 2 , x 3 ) (A. 6) 

Since x 3 on the surface is not an independent variable, 
therefore , 

g = qtx^, x 2 , x 3 (x ir x 2 )] = qtx.^ x 2 ) (A. 7) 

The relation between the partial derivatives of the flow property q 
in the three-dimensional field with those on the surface can 
be written as: 


n„ 


i3_ = 13 I 


ax. 


3X, 


n 3 h 3 


1 ag_ 

3x- 


3q , 3q ^ n 2 h 2 3g (A 

3x 2 3x 2 n 3 h 3 8 x 3 

where n^, n 2 , and n 3 are the components of a unit vector n that 
is normal to the prescribed surface in the x^, x 7 and x 3 
directions , respectively . 
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The previous derivation is valid for any generalized curvi- 
linear coordinate system. Following Wu [ 1] , the best coordinate 
choice for the surface are the meridional distance m, the 
tangential angle and the normal distance n to the surface 
(see Fig. 1) . 

The following relations apply for the (m-$-n) coordinate 
system 


n l “ n 2 = 0 ' 



(A. 9) 


and the metric or Lamie coefficients h^, h^, h^ are given by: 


h^ = 1 , h 2 = r , h^ = 1 (A. 10) 

In the present study, the number of the surfaces, to which the 
passage is divided, is chosen to be large. Therefore, the 
filament thickness b or each surface is considered to be small 
compared to the radius r. Consequently, for those surfaces which 
are located away from the hub and shroud boundaries (Fig. 3) , 
one can consider that the change of flow properties across the 
filament thickness, b, is neglected. Thus 


3q _ 3q 
3^3 3n 3 


0 


(A. 11) 


Using equations (A. 9), (A. 10) and (A. 11) to rewrite the right 

hand side of equation (A. 9), we obtain 


3q 3q j[C£ 

3x, 3x, 3m 


3q _ j3£ 
3x 2 3<f> 


ax. 


R^ESDUCffilLlK^raE 
OKICtW- is VQm 


(A. 12) 


(A. 13) 


Substituting the expressions in equations (A. 11), (A. 12) and (A. 13) 
for the derivatives of all quantities in the equations of motion 
(A.l), (A. 2) and (A. 3), we obtain: 
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Conservation of Mass 


If + Is + ? k <P V + ? M m H = 0 


(A. 14) 


where W and W. are the components of relative velocity vector 
III 9 

in the meridional and tangential direction respectively. 
Conservation of Momentum in the Meridional Direction 


3W 3W W. 3W w 2 - 

» rr + + rsina - 2 !i V lm) 


= „ ^ - 1 iL_ r u ( 


3 <fj r 

aw w 


m 


3m 3 3m e 3 i 


i 3W* 

- sxna -r - y) ] 


, . aw a aw. w, , aw m 

-j* - {§— ( 2 y r 5 —^) + 7 r-r-[y (r— — sina + — 3 

r 3m e 3m 3<p e 3m r r d<p 


u a , aw. w m 

2 -S. [i — 4 . + -SL sina] sina 

r l r 3$ r 


(A. 15) 


where a is the angle between the meridional streamline and the 

dr 

axis of rotation (z) , as shown in Fig. 1, such that sina = 
Conservation of Momentum in the Tangential Direction 


3W 3W 3W 

pr + P (rW t— + W, + W W sina + 2JIr W_ sina) 

3 1 m am 9 o<p m 9 iu 

, a aw^ w M , aw, 

3$“3 34 ^e 3m ^ r n + r 34> }3 

- aw, aw - . aw, w 

+ 3S tu e (E Sr- W * sina + 3r )I * 3* t2u e ( S^- + F !1 sinc0 1 


3W W 

+ “e sina [ a 5 T “ V- sina + I I? V 


(A. IS) 
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Adopting equations (A. 14) , {A. 15) and (A. 16). leads to a system 
that contains derivatives with respect to two space variables 
only, namely m and Consequently, the flow over the surface 
may be treated in a two dimensional manner. 

In order to obtain an equation for a stream function, 
reference [1] introduces the concept of an integrating factor, b, 
such that the equation of continuity becomes: 


br H + k (br ■> V + 1* < b » V - 0 


(A. 17) 


where the factor, b, satisfies the following relation. 


W. , 

m Sb , _jp 9 t- l - n 
b 3m br 3<p “ u 


(A. 18) 


A stream function may now be defined such that 


_ M 1 
W m ~ b pr 


3 <j> 



M 1 3* 
b p 3m 


(A. 19) 


The previous equations indicate that the integrating factor, b, 

is nothing more than the filament thickness of the surface S^. 

It is worth noting that the above system of equations are 

valid for any computational surface providing that its geometry 

has been defined. The system can be applied in the rotating frame 

as stated, or in the absolute stationary frame by setting a = 0. 
ctx 

When ^ = 1 or o = 90°, the equations represent the flow in a pure 
radial machine. Also, when ^ = 0 or a = 0°, the equations 
represent the flow in an axial machine,’’ 


REPRODUCIBILITY of THE 1 

CRlGpv T AL PAGE IS POOR^ 


80 


Energy Equation 

For steady relative flow through a turbomachine rotor, the 
conservation of energy given as equation (3) in Section 1, can 
be rewritten in the following form: 

p e 

pW* Vh = W*Vp + D + V (~ Vh) (A. 20) 

where Pr is the Prandtl number and D is the energy dissipation 
function. 

It is more convenient to express the energy equation {A. 20) 
in terms of the total enthalpy H of the gas, and its velocity 
components- The total enthalpy for turbulent flow is 

expressed as follows. 

H = h + W 2 /2 + OW^r + fl 2 r 2 /2 + E (A. 21) 

where E is the kinetic energy of turbulence. 

Using the above equation, the energy equation (A. 20) reduces 

to: 

2 2 2 

pW*VH = pW-V[“4* flW^r + ~ 2 r ■ + E] + W*Vp 

+ v * [pf V(H-| S )J + D (A. 22) 


Equation (A. 22) may be expressed in a slightly different form by 
eliminating the pressure term in the right hand side using the 
momentum equation (1) of Section 1, with the following result 

pW* VH = pW*7lQW^r + n 2 r 2 4- E] 

u e W 2 o 2 ^ 2 

+ ’< h - 2- * n V - V- - E >I 

- W*V (u w) + iw*V(jj V«W) + D (A. 23) 

S j 0 
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It is to be noted that in formulating the above relation, use 

has been made of the following vector identities: 

____ w 2 

W-flxw = 0 , w*w*v w = w-v ~ 

2 

W-R = W-V ~~ , VxW = a 


where a is the vorticity vector. 

In order to preserve simplicity of concept and for better 
organization, it would be very advantageous, as pointed out in ref. 
[22] , to replace the kinetic energy of turbulence convective term, 
pW-VE, that appear in the above equation by its equivalent diffusion 
expression. The exact form can be expressed, using the transport 
equation of the kinetic energy of turbulence, as follows: 


pW'VE 



VE) 


(A. 24) 


where S CE has the significance of Schmidt number for the kinetic 
energy of turbulence. It is worth noting that the adoption of 
equation (A. 24) implies that the generation and decay terms in the 
full kinetic energy of turbulence equation are in balance. 

Substituting now equation (A. 24) in (A. 23) , we obtain: 


J-i 2 

pW-VH = pw-v[nw,r + fl 2 r 2 ] + V-[-^ V(H - 

9 Pit 2 <p 

$ 2 2 r 2 - 

2 E) + s~“ VE] ~ W-7x(u u) 

CE e 

4 - _ 

-J- y W-V(U e V*W) + D 


(A. 25) 


For the (m,<ji,n) coordinate system, associated with the computa- 
tional surface of Fig. 1, the 7 operator can be expressed as 


V 


- a e $ a 

e m Bm r 3<j> 


(A. 26) 


where e m and e, denote unit vectors in the meridional and 
m (j> 

tangential directions respectively. The velocity vector W, and 
the vorticity vector to, may be written as follows in this 
system of coordinates: 


W = We + 
m m 


V* 


_ M 1 ,1 3* - 
bp l r 3$ e m 


ii i ] 

3m <f> J 


to = ton 


(A. 27) 


When the relations (A. 26) and {A. 27) are used to fully expand the 
energy equation (A. 25), the following result is obtained. 


b C 3m (H r ||)- i L.& 






+ {y r [— - w 12 . _/i 1 s 3E. , 

3m lM e l Pr 3m 


S cs Pr' 3 nr* 


+ |x {-& [i- jgC/2. _ / 1 1 , 3E n rT 3 

3< f> r L Pr 3d> ( s^ 5?) axl > - W A 2_ 


3*'" 'S CE p r J 3$^ “ W <J, <V a) 


W . 
m 3 

r 3$ 


(urn) - br + G- - 0 


(A. 28) 


where 

G. 


- - a | HsHV* 20 ' §* [< V + ^> £» 


3U)i 


3 r K e _ 3 ^ 1+ 1 i_rl^ L_ (w r+ 2§1) ]} 

? J 1 r 34> L Pr 3$ * 2 


+ [— r ~ (W.r+ 

3m L Pr 3m <£ 


(A. 29) 
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and the dissipation, function , 


3W 2 , 3W, 

D = 2u e * ( 3m“ } + a$~” 


+ U e C 3m 


3W, , 3W 

g 4 _ 1 in 

r 3g 


D, is given by 

W ~ 

+ — sina) “} 
r 

W <f> . \ 2 

- ™ sma} 


(A. 30) 
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APPENDIX B 


DERIVATIVE TRANSFORMATION 


This appendix presents the mathematical expressions used to 
transform the partial derivatives from the physical space of 
(m,4i) coordinates to the computational domain of (£,n) coor- 
dinates. A few relations involving directional derivatives 
either normal or tangent to a line of constant £ as well as a 
line of constant n are also included. Since the purpose of 
this appendix is to provide a quick reference only, most of 
the algebraic development is omitted. 


Partial Derivative Transformation 

As pointed out earlier in Section 2, two transformations are 
employed to transform the physical space of Fig. 5a into the 
unit square of Fig. 5c. These transformations are given by 
equations {31} and (35) . The partial derivatives of any scalar 
function, f, are transformed utilizing the chain rule as 
follows: 


3f _ 5f 3x f 3f 3<ft _ 1 3 f 
3m '3x 3m 3<j> 3m r 3x 


1 1 li + M £lh 

r ^3£ 3x an 3x J 

1 ,34> 3f d± 3f. 
rJ l 3n 35 ~ 35 3n J 


(B.l) 


3f 3 f 9 4> , 3f 3n 

3(f> 3£ 3£ + 3n 3<t> 

_ , 3x _3j: 3x 3_f. /T 

l 3£ 3n " 3n 35 J/ (B . 2) 


85 


where J is the Jacobian of transformation given by: 

-r ... 3x 3 6 _ 3x 3<j> (b. 3) 

u 34 3 n 3 n 34 

Likewise , higher order derivatives, which appear in the flow 
governing equations, are transformed using the following 
relations : 


r L. £r M) = iff 

r 9m lr 3m J 2 


= i__ rfii^ 2 iff _ o 34* 3* 9 2 f . ,36*2 3 2 f, 

? a? 2 2 H In a| 5 n * ( 3 ? } ^ 2 ] 


* L_ rr^ 2 l 2 ! _ , 3* 34* a 2 <f> , ,36*2 3 2 6-, 

J 3 L ^n J 3^2 2 ST air + { W ^ 


r 9x 9f 3x 3f * t r ,36* 2 3 2 x n 36 36 3 2 x 

»*• ir Tq'J + l T* “ ^ 


•9n 3? 34 3n 


l 3 n 


ar 


+ rli) 2 if*] fit 3f . 14 3f * •> 
^34 J g n 2 Jt 34 3n an 34 )} 


34 3n 343n 


(B.4) 


a 2 f _ 1 r ,3x 2 3 2 f o 3x 3x 3 2 f 

n 2 n ^ r» * - 


36 


3 C J 


34 3n scan 

2 


+ (!*) 2 ifl] 

‘ 3 5’ a, 23 


+ i_ rrf 3 x i 2 ^ - 3x 3x • ' ax ' 2 a**, 

+ J 3 U( 3 lT ) ^2 2 a? 3 ? 343 TT + ( 34 } ^ 2 ] 


fii£ ii _ 3x 3_f . r .3x 2 3 2 x _ , 3x 3x 3 2 x 

3n 34 34 3n J n 3n J * r 2 2 34 3n 343n 


, 3x* 2 j3_x- .36 3f 36 3f 

34 J g n 2 J 3 4 3n 3n 34 


34‘ 
■) 


(E. 5) 


36 


- — (r — 
3<i> ir 3m ; 


2 

3 £ = i_ r . ax a<&. a 2 f 

3(|)3x 2 U 35 3 ti 3n 35' 3£3n 


3x 36 3 f 

3? 35 an 2 


ax 36 a' 
an an . 


as 


I] + I f 

2 J + 3n 


r l f 3x 3t*> 3 J 

L 3 

J 


-— J^r. — _ _ 3x a_6 3J\ 

a? a? an an 35 35 


M 1 ^ _ 9 _X 3 £ 

J 2 911 7? 


ax a 2 6 


35 353r| 


•)] + 


3f r l 
35 


/ ax a 6 a j 
l an an 35 


3x 36 3J. 

35 an an' 


7 


,3x 3 $ 

35 87 ' 


ax ‘a 2 ® 


an 35an 


)] 


(B.6) 


Directional Derivatives 7 

In most fluid flow analyses, the directional derivatives of 
the flow properties, along or normal to the boundaries, are 
often needed to evaluate various boundary conditions. These 
quantities can be easily obtained if the unit vectors tangent 
or normal to a line of constant ^ or n, that coincide with the 
physical boundaries, are specified. As pointed out in refer- 
ence [12] , these vectors can be expressed in terms of the 
transformation parameters f^, 14, etc. as follows: 

, d c, d H d £ 

!e t' 5 S m + 3? e <t> ) ^ 

l h\ = ( !f + H V /,/? 


t Some of the relations given in this appendix can be found in 
reference [12] . They are repeated here howver for completeness 
and easy reference for the development of the governing 
equations . 
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(fi—) £• 
n £ 


<IK ' W %>/* 


{e ) 

n n 


<" It 5 m + If 


(B.7) 


Where (e.) - and (e ) represent the unit vectors tangent to 
^ r| ^ 

constant £ and constant rt lines respectively. While, (e^) ^ 

and (e ) denote the unit vectors normal to the constant £ and 
n n * 

the constant n lines . The above mentioned unit vectors are 
shown in Fig. B-l as they appear in the physical domain. 

The transformation parameters 6 and y in equation (B.7) 
are given by 


6 - ,S± 2 

6 " ( 3n } ( 3n 


Y = (M ) 2 + { li } 2 


The directional derivative of any scalar function, f, in any 
direction e is given by 


TC - < 5 • 7)£ 


where Vf can be written as 


_ L_ r (ii if _ li M\ S + ^ ^ l£] «1 

~ rJ u 3n 3£ 3£ 3n ; m l 3£ 3n 3n 3£ J e <p J 


(B.8) 


Associating the direction e with the unit vectors normal and 
tangent to lines of constant £ or rir we have 


7^: = He ) r *V]f = 

8(e t ) 5 t; 


n 5 t V*« - «//r 
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If - 8 ff>A"v 


where 

o — 3 4> 94* 

B 3? 3n 35 3n 


(B . 9 ) 


(B. 10) 


With these relations th s different derivative boundary conditions 
associated with the flow governing equations may be evaluated 
in a straightforward manner, providing that the values of the 
different transformation parameters are available. 
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